BLOCK KRYLOV SUBSPACE EXACT TIME INTEGRATION 
OF LINEAR ODE SYSTEMS 
PART 1: ALGORITHM DESCRIPTION* 
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Abstract. We propose a time-exact Krylov-subspace-based method for solving linear ODE (or- 
dinary differential equation) systems of the form y' = —Ay+g{t), where y{t) is the unknown function. 
The method consists of two stages. The first stage is an accurate polynomial approximation of the 
source term g(t), constructed with the help of the truncated SVD (singular value decomposition). 
The second stage is a special residual-based block Krylov subspace method. 

The accuracy of the method is only restricted by the accuracy of the polynomial approximation 
and by the error of the block Krylov process. Since both errors can, in principle, be made arbitrarily 
small, this yields, at some costs, a time-exact method. 
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1. Problem formulation. Consider initial- value problem (IVP) 

U/ = -Ay + g{t), 
\y{0)^v, t e [0,T], 

where y{t) is the unknown vector function, y : M — > M", and the matrix A G R"^", 
vector function g : R — > M", and vector v € E" are given. 

Let y{t) = y{t) — v (meaning that y{t) — y{t) —v for all t). Note that the function 
y{t) satisfies IVP 

[y' = -Ay + m, .^2) 
\y(0) = 0, te [0,T], 

where g{t) = g{t) — Av. We will assume that the IVP is brought to the equivalent 
form (jl.2p and, for simplicity, we omit the tilde sign ~ in l|1.2p . 

2. Polynomial approximation. We now describe the first stage of the method, 
the best fit polynomial approximation of the source term g{t). Choose s points 
= ti < <2 < • ■ ■ < ts-i < ts = T on the time interval [0,T]. The polynomial 
approximation is based on the truncated SVD (singular value decomposition) of the 
matrix 

G^[gih) gih) ... g(^,)]eR"^^ 

whose columns are samples g{ti), i — l,...s, of the vector function g{t). More 
precisely, let 

G^UtV^, S = diag(ai,...,a,) GR^^^ di ^ . . . ^ a, ^ 0, (2.1) 
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be the thin SVD [3 Section 2.5.4], where the matrices U e M"^" and V G R"""" have 
orthonormal columns Ui, . . . , Ug and vi, . . . , respectively. An approximation to 
G can be obtained by truncating the SVD as 

s m 

G = iltV^ = ^ a^u.vj w ^ cTjUjwf = UYy^, m< s, (2.2) 

i=l 4=1 

where S e K™^™ = diag(cri, . . . , ct™) and the matrices U e R"^™ and V G R"""" 
are formed by the first m columns of U and V, respectively. Denote the obtained 
approximate matrix by G = UT,V^ . If follows from (|2.2I) that the SVD of G — G is 
readily available as '^i^m+i '^i'^i'^I' ■ Hence, for the 2-norm and Frobenius norm of 
the error G ~ G holds [31 Section 2.5.3]: 

||G - GII2 = cr„+i, IjG - G||^ = cr^^i H hfTs- 

Looking at SVD identity (12. ip columnwise, we see that every sample value g{ti) of 
the function g(t) can be approximated by a linear combination of the vectors ui, . . . , 

g{ti) = {(TiVii)ui + {(T2Vi2)u2 H h {asVis)Us 

W {aiVii)ui + {(T2Vi2)u2 H h icrsVim)Um, 

where Vij are the entries of the unitary matrix V. Following the approach of [2], 
we consider the coefficients of these linear combinations, namely cTjVij, j — 1, . . . , m, 
as values of some unknown functions fj{t) at ti. These functions can be easily ap- 
proximated, at a low cost (typically m n) and with a very high accuracy, by a 
polynomial fit [2]. This yields 

9{U) ~ fl{U)ui + f2(ti)u2 H h Jm{ti)Um 

K. pi{ti)ui +P2{ti)u2 H \-Pmiti)Um- 

For simplicity, we assume that all the best-fit polynomials have the same the order 
r. Packing the polynomials Pj{t), j = 1, . . . ,m, in one polynomial vector function 
p{t) = {pi{t), . . . ,pm{t))'^ , we obtain a polynomial approximation 



g{t) « Up{t). (2.4) 

There are three sources contributing to the approximation error here. First, the 
quality of the approximation is influenced by the choice of the sample points , . . . , 
ts- Second, by the number of terms m in the SVD truncation ()2.2|) and, finally, by 
the polynomial best fit in (|2.3p . All these errors can be easily controlled when the 
approximation is constructed [2 , thus giving possibility for an adaptive approximation 
procedure. With (|2.4I) . the original initial- value problem (|1.2p takes the form 



(y'^ -Ay + Up{t), 
\y(0)=0, ie[0,T], 

We now introduce a block Krylov subspace method to solve this problem. 
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3. Residual-based block Krylov subspace method. To construct a Krylov 
subspace block iterative method for solving (|2.5|) . we use the exponential residual 
concept described in [T]. Choosing the initial guess yo{t) to be a zero vector function, 
we see that the corresponding initial residual is 

ro(i) = -Ayoit) - y^(i) + C/p(i) = C/p(t). (3.1) 

We follow the approach of [1 , where the approximate solution ykit) at Krylov iteration 
k is obtained as 

Here the vector function ^fe(t) is the Krylov subspace approximate solution of the 
correction problem 

\^(0) = 0, [0,T], 

Note that if ^fc(t) solves p.2p exactly then yk{t) is the sought-after exact solution 
of (|2.5p . We solve (|3.2p by projecting it onto a block Krylov subspace defined as 

ICkiA, U) EE span [U, AU, A^U, A''~^U} , 

with dimension at most k-m. An orthonormal basis for this subspace can be generated 
by the block Arnoldi or Lanczos process described e.g. in O |4] . The process produces, 
after k block steps, matrices 

Here e M"""", Vi is taken to be the matrix U produced by the truncated SVD (|2?2)) 
and V[fe+i] has orthonormal columns spanning the Krylov subspace, namely, 

colspan( V[j.] ) = Kk(A, U). 

The matrix H\j.^i^k] is block upper Hessenberg, with mxm blocks Hij, i — 1, . . . , fc+1, 
j = 1, . . . ,k. The matrices V[/j+i] and H[k+i,k] satisfy the block Arnoldi (Lanczos) 
decomposition [51 [J 

AV[k] = V[k+i]H[k+i,k] = V[k]H[k.k] + Vk+iHk+i,k^k ^ (3-3) 

where Hk+i^k is the only nonzero block in the last k + 1 block row of H^k+i.k] and 
Ek S E"^*^ is formed by the last m columns of the km x km identity matrix. 
The Krylov subspace solution (t) is computed as 

6(t) = V{k]u{t), 
where u{t) solves the projected IVP 

(u'{t) = -H[k,k]u{t) + V[l]roit), 
\u(0) = 0, te[0,T]. 

Note that 
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where Ei G ^kmxm formed by the first m columns of the km x km identity matrix. 
Using (I3.ip . (I3.3P and p.4p . we see that for the exponential residual rk{t) of the 
solution yk{t) holds 

rfc(i) = -Ayk -y'k + Up{t) = -Ayo -y'o~ AVik]u{t) ~ V[fc]w'(t) + Up{t) = 
= ro{t)-AV[k\u{t)-Vyk]u'{t) = 

= ro{t) - {V[k]H[k,k] + Vk+iHk+i,kEl)u{t) - V[k]u'{t) = (3.5) 

= ro(t) - V[fe](i/[fe,fe]ii(t) - Vk+iHk+i^kElu{t) = 

= ro(i) - T4£;ip(t) - Vk+iHk+i,kElu{t) = -Vk+iHk+i,kElu{t). 

A similar expression for the exponential residual is obtained in [P for a non-block 
Krylov subspace method. There are two important messages relation p.5|) provides. 
First, the residual can be computed efficiently during the iteration process because 
the matrices Vk+i and Hk+i^k are readily available in the Arnoldi or Lanczos process. 
Second, the residual after k block steps has the same form as the initial residual (|3.ip . 
namely it is a matrix of m orthonormal columns times a time dependent vector func- 
tion. This allows for a restart in the block Krylov method: set yoit) := yk{t), then 
relation p.ip holds with U := Vk+i and p{t) :— —Hk+i,kE'^u{t). The just described 
correction with k block Krylov iterations can then be repeated, which results in a 
restarted block Krylov subspace method for solving (|2.5p . 
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